function SS=SteadyState_Solver(PP,PSI_0,DiscFactor_0)
% clear;
% clc;
% PP          =   Setup_PP(struct());
%% Preliminaries
if PP.InfoPrint>3
    PrintFlag   =   1;
else 
    PrintFlag   =   0;
end
if nargin<=1
    PSI_0       =   1.2;
    DiscFactor_0=   0.1;
end

%% Initiate the SS
if PrintFlag
    fprintf(['Start Solving the Steady State ...\n']);
end
TimeStart   =   tic;
SS          =   struct();
%--------------------------------------------------------------------------
% Exogenously Fixed
%   Note: 
%--------------------------------------------------------------------------
SS.Pir      =   PP.Pir_SS;

SS.N        =   PP.N_SS;
SS.wN       =   1;
SS.ir       =   PP.i_SS;
SS.rr       =   SS.ir-SS.Pir;

SS.TAU      =   PP.TAU_SS;
SS.T        =   PP.T_SS;

%--------------------------------------------------------------------------
% Approximation Structure for Value Function, Policy Function, and Distribution
%   Note: 
%--------------------------------------------------------------------------
[SS.FunApp,SS.PolApp] ...
            =   SteadyState_FunPolApp(PP,SS);
SS.DistApp  =   SteadyState_DistApp(PP,SS);


%% Solve the SS by Finding to Equilibrium 
%--------------------------------------------------------------------------
% Solve the root by Newton Update 
%   Note: 
%       1. Save time by recycling the UCoef from the intermediate step.
%       2. Fastest, but might not be robust when Diff(PSI,ra) is too flat.
%--------------------------------------------------------------------------
BellmanItSetup=   struct('ValItNum',0,'ColItNum',1,'TolItNum',50,...
                         'PrintFlag',PrintFlag,'DampenCoef',1);
Input_0     =   [PSI_0;DiscFactor_0];
[SS,ExitFlag]=  SteadyState_Solver_Newton(PP,SS,BellmanItSetup,Input_0);
%--------------------------------------------------------------------------
% Recollect the Steady State
%   Note: 
%--------------------------------------------------------------------------
TimeElapsed     =   toc(TimeStart);
if PrintFlag
    fprintf(['Steady State is Solved, Elapse Time=',num2str(TimeElapsed,'%.2f'),'s','\n',...
             'Equilibrium Residual: \n',...
             'Labor Market, Total Labor Quantity: ',num2str(SS.RES(2),'%.2e'),', \n', ...
             'Domestic Bond Market: ',num2str(SS.RES(1),'%.2e'),', \n',...
             'Final Good Market: ', num2str(SS.RES_FinalGood),'.\n']);
end

%% State Space Reduction
%--------------------------------------------------------------------------
% Distribution
%--------------------------------------------------------------------------
[SS.Dist,SS.CopulaF] ...
            =   DimReduction_Distribution(PP,SS,'Zip',SS.Dist_fz.QW_fz);
%--------------------------------------------------------------------------
% Value Function
%--------------------------------------------------------------------------
SS.ValFun   =   SS.UCoef-SS.UCoef;

%% Important Statistics
SS.AggStat  =   [SS.Agg.Avg.c.poor;SS.Agg.Avg.c.rich];

%% Residuals

